grb_batse = read.csv("grb_dataset.csv")
head(grb_batse)
grb_dur = read.csv("BATSE-4Br-dat.csv")

head(grb_batse)
head(grb_dur)

summary(grb_batse)
summary(grb_dur)

grb <- merge(grb_batse, grb_dur, by = "TrigNo", all = TRUE)
head(grb)
grb = grb[,-c(1:8,11)]
head(grb)

library(circular)
library(CircStats)

watson.test(grb$GLON, alpha = 0.01, dist = "vonmises")
watson.test(grb$GLAT, alpha = 0.01, dist = "vonmises")

library(skmeans)
library(Matrix)
library(slam)
dim(grb)
grb = na.omit(grb)

x = as.simple_triplet_matrix(grb)


library(skmeans)
library(slam)

# Assuming x_clean is your cleaned simple triplet matrix
k_range <- 2:10
obj_values <- numeric(length(k_range))

# Compute objective values
for (i in seq_along(k_range)) {
  k <- k_range[i]
  skm <- skmeans(x, k = k)
  obj_values[i] <- skm$value
}

# Compute first differences (slope)
first_diff <- diff(obj_values)

# Compute second differences (curvature)
second_diff <- diff(first_diff)

# Find the k with the maximum curvature (largest drop in improvement)
# We look at the absolute value of second_diff since curvature could be negative
k_elbow <- k_range[which.max(abs(second_diff)) + 1]  # +1 because second_diff is shorter

# Plot elbow curve with optimal k marked
plot(k_range, obj_values, type = "b", pch = 19, 
     xlab = "Number of clusters (k)", ylab = "Within-Cluster Cosine Similarity",
     main = "Elbow Method with Optimal k")
lines(k_range, obj_values, col = "blue")
abline(v = k_elbow, col = "red", lty = 2)
cat("Optimal number of clusters (max curvature):", k_elbow, "\n")

# Optional: Plot second differences to confirm
plot(k_range[2:(length(k_range)-1)], abs(second_diff), type = "b", pch = 19,
     xlab = "Number of clusters (k)", ylab = "Absolute Second Difference",
     main = "Curvature Plot")
abline(v = k_elbow, col = "red", lty = 2)

### As per this we get that 3 is the optimal number of clusters
set.seed(1423)
hpart = skmeans(x, 3)
hpart$cluster

grb$cluster = hpart$cluster
head(grb)


# Load required libraries
library(MASS)  # For LDA
library(ggplot2)  # For visualization

# Assuming grb is your dataset
# Select features used in skmeans (adjust columns as needed)
features <- grb[, c("GLAT","GLON", "T50", "T90", "F1", "F2", "F3", "F4", "P64", "P256", "P1024")]
clusters <- grb$cluster


# Normalize rows to unit length (spherical assumption)
row_norms <- sqrt(rowSums(features^2))
features_norm <- features / row_norms
# Check the first few rows
head(features_norm)

# Check variance within each cluster
for (i in 1:ncol(features_norm)) {
  var_by_cluster <- tapply(features_norm[, i], clusters, var)
  cat("Variable", colnames(features_norm)[i], "variance by cluster:\n")
  print(var_by_cluster)
  if (any(var_by_cluster == 0, na.rm = TRUE)) {
    cat("  -> Zero variance detected in cluster(s):", which(var_by_cluster == 0), "\n")
  }
}

# Load libraries
library(MASS)
library(ggplot2)

# Prepare data (assuming grb includes GLON and GLAT as features)
features <- grb[, c("GLON", "GLAT", "T50", "T90", "P64", "P256", "P1024")]
clusters <- grb$cluster

# Normalize rows to unit length
row_norms <- sqrt(rowSums(features^2))
features_norm <- features / row_norms

# Check for any remaining issues
print(head(features_norm))
print(table(clusters))  # Ensure no clusters have < 2 points

# Perform LDA
lda_model <- lda(features_norm, grouping = clusters)
print(lda_model)

# Predict and evaluate
lda_pred <- predict(lda_model, features_norm)
table_true_pred <- table(True = clusters, Predicted = lda_pred$class)
print(table_true_pred)
accuracy <- sum(diag(table_true_pred)) / sum(table_true_pred)
cat("Classification accuracy:", accuracy, "\n")

# Visualize
lda_scores <- as.data.frame(lda_pred$x)
lda_scores$Cluster <- as.factor(clusters)
ggplot(lda_scores, aes(x = LD1, y = LD2, color = Cluster)) +
  geom_point(size = 3) +
  labs(title = "LDA Projection of Spherical Clusters (Reduced Features)",
       x = "LD1", y = "LD2") +
  theme_minimal()






#### Spherical Model Based Clustering


# Function to compute the vMF density (log-scale for stability)
vmf_density <- function(x, mu, kappa) {
  d <- ncol(x)  # Dimensionality
  log_C <- (kappa^(d/2 - 1)) / ((2*pi)^(d/2) * besselI(kappa, d/2 - 1))  # Normalization constant
  log_density <- log_C + kappa * (x %*% mu)  # x and mu are unit vectors
  return(log_density)
}

# EM algorithm for vMF mixture
vmf_mixture_em <- function(X, k, max_iter = 100, tol = 1e-6) {
  n <- nrow(X)  # Number of observations
  d <- ncol(X)  # Dimensionality
  
  # Initialize parameters
  set.seed(123)
  pi <- rep(1/k, k)  # Equal mixing proportions
  mu <- t(apply(matrix(runif(k * d, -1, 1), k, d), 1, function(x) x / sqrt(sum(x^2))))  # Random unit vectors
  kappa <- rep(10, k)  # Initial concentration (can tune this)
  z <- matrix(0, n, k)  # Responsibilities
  
  log_lik_old <- -Inf
  
  for (iter in 1:max_iter) {
    # E-step: Compute responsibilities
    for (j in 1:k) {
      z[, j] <- log(pi[j]) + vmf_density(X, mu[j, ], kappa[j])
    }
    z <- exp(z - apply(z, 1, function(row) max(row)))  # Subtract max for stability
    z <- z / rowSums(z)  # Normalize to probabilities
    
    # M-step: Update parameters
    # Mixing proportions
    pi <- colMeans(z)
    
    # Mean directions (normalized weighted average)
    for (j in 1:k) {
      r_j <- colSums(z[, j] * X)  # Weighted sum of observations
      mu[j, ] <- r_j / sqrt(sum(r_j^2))  # Normalize to unit vector
      # Update kappa (approximation via maximum likelihood)
      r_bar <- sqrt(sum(r_j^2)) / sum(z[, j])  # Average resultant length
      kappa[j] <- (r_bar * d - r_bar^3) / (1 - r_bar^2)  # Approximation for kappa
    }
    
    # Compute log-likelihood
    log_lik <- 0
    for (i in 1:n) {
      ll_i <- numeric(k)
      for (j in 1:k) {
        ll_i[j] <- log(pi[j]) + vmf_density(X[i, , drop = FALSE], mu[j, ], kappa[j])
      }
      log_lik <- log_lik + log(sum(exp(ll_i - max(ll_i)))) + max(ll_i)  # Log-sum-exp trick
    }
    
    # Check convergence
    if (abs(log_lik - log_lik_old) < tol) break
    log_lik_old <- log_lik
  }
  
  # Cluster assignments
  clusters <- apply(z, 1, which.max)
  
  return(list(pi = pi, mu = mu, kappa = kappa, z = z, clusters = clusters, log_lik = log_lik))
}
dim(grb)


X <- as.matrix(grb[,-12] / sqrt(rowSums(grb[,-12]^2)))  # Normalize to unit sphere
dd = X

# Run the EM algorithm with 3 clusters
k <- 3
set.seed(12433)
result <- vmf_mixture_em(X, k)
result
grb$vmf_cluster = result$clusters

# Results
cat("Mixing proportions:", result$pi, "\n")
cat("Concentration parameters:", result$kappa, "\n")
cat("Log-likelihood:", result$log_lik, "\n")

# Perform LDA
lda_model_modelB <- lda(features_norm, grouping = result$clusters)
print(lda_model_modelB)

# Predict and evaluate
lda_pred <- predict(lda_model_modelB, features_norm)
table_true_pred <- table(True = result$clusters, Predicted = lda_pred$class)
print(table_true_pred)
accuracy <- sum(diag(table_true_pred)) / sum(table_true_pred)
cat("Classification accuracy:", accuracy, "\n")


# Visualize
lda_scores <- as.data.frame(lda_pred$x)
lda_scores$Cluster <- as.factor(result$clusters)
ggplot(lda_scores, aes(x = LD1, y = LD2, color = Cluster)) +
  geom_point(size = 3) +
  labs(title = "LDA Projection of Spherical Model Based Clusters (Reduced Features)",
       x = "LD1", y = "LD2") +
  theme_minimal()

library(mclust)
skmeans_clusters = hpart$cluster
vmf_clusters = result$clusters

cont_table = table(skmeansss = skmeans_clusters, vMFc = vmf_clusters)
cont_table

ari = adjustedRandIndex(skmeans_clusters, vmf_clusters)
ari

# Assign T90-based classes
grb$T90_class <- ifelse(grb$T90 < 2, "Short",
                        ifelse(grb$T90 <= 10, "Intermediate", "Long"))

# Create comparison tables for Spherical K-means and vMF
skmeans_t90_table <- table(Cluster = skmeans_clusters, T90_Class = grb$T90_class)
vmf_t90_table <- table(Cluster = vmf_clusters, T90_Class = grb$T90_class)

# Print tables
print("Spherical K-means vs. T90 Classes:")
print(skmeans_t90_table)
print("vMF vs. T90 Classes:")
print(vmf_t90_table)


# Compare with Horvath (1998b) classifications (simulated here; replace with actual data if available)
# For demonstration, assume Horvath's classes are similar to T90-based
horvath_classes <- grb$T90_class  # Placeholder; replace with actual Horvath (1998b) labels
skmeans_horvath_agreement <- mean(skmeans_clusters == as.numeric(factor(horvath_classes)))
vmf_horvath_agreement <- mean(vmf_clusters == as.numeric(factor(horvath_classes)))
cat("Agreement with Horvath (1998b) - Spherical K-means:", skmeans_horvath_agreement, "\n")
cat("Agreement with Horvath (1998b) - vMF:", vmf_horvath_agreement, "\n")





####### Kmeans Clustering #####
set.seed(122222)
head(grb)
km = kmeans(grb[,-c(12:14)], 3)
# Perform LDA
grb$km_cluster = km$cluster

lda_model_km <- lda(features_norm, grouping = km$cluster)
print(lda_model_km)

# Predict and evaluate
lda_pred <- predict(lda_model_km, features_norm)
table_true_pred <- table(True = km$cluster, Predicted = lda_pred$class)
print(table_true_pred)
accuracy <- sum(diag(table_true_pred)) / sum(table_true_pred)
cat("Classification accuracy:", accuracy, "\n")

# Visualize
lda_scores <- as.data.frame(lda_pred$x)
lda_scores$Cluster <- as.factor(km$cluster)
ggplot(lda_scores, aes(x = LD1, y = LD2, color = Cluster)) +
  geom_point(size = 3) +
  labs(title = "LDA Projection of kmeans Clusters (Reduced Features)",
       x = "LD1", y = "LD2") +
  theme_minimal()

##### GMM ######
library(mclust)
set.seed(12222231)
gmm = Mclust(grb[,-12], 3)



# Create comparison tables for Spherical K-means and vMF
kmeans_t90_table <- table(Cluster = km$cluster, T90_Class = grb$T90_class)
gmm_t90_table <- table(Cluster = gmm$classification, T90_Class = grb$T90_class)

print(kmeans_t90_table)

print(gmm_t90_table)



kmeans_horvath_agreement <- mean(km$cluster== as.numeric(factor(horvath_classes)))
gmm_horvath_agreement <- mean(gmm$classification == as.numeric(factor(horvath_classes)))
cat("Agreement with Horvath (1998b) -  K-means:", kmeans_horvath_agreement, "\n")
cat("Agreement with Horvath (1998b) - gMM:", gmm_horvath_agreement, "\n")


lda_model_mbc <- lda(features_norm, grouping = gmm$classification)
print(lda_model_mbc)

# Predict and evaluate
lda_pred <- predict(lda_model_mbc, features_norm)
table_true_pred <- table(True = gmm$classification, Predicted = lda_pred$class)
print(table_true_pred)
accuracy <- sum(diag(table_true_pred)) / sum(table_true_pred)
cat("Classification accuracy:", accuracy, "\n")

# Visualize
lda_scores <- as.data.frame(lda_pred$x)
lda_scores$Cluster <- as.factor(gmm$classification)
ggplot(lda_scores, aes(x = LD1, y = LD2, color = Cluster)) +
  geom_point(size = 3) +
  labs(title = "LDA Projection of Mixture Model Based Clusters (Reduced Features)",
       x = "LD1", y = "LD2") +
  theme_minimal()










#### Plot over a sphere ####
library(rgl)

clusters = hpart$cluster
# Step 4: Plot the results on a 3D sphere
# Open a 3D plotting window
open3d()
library(Directional)

ss = data.frame(grb$GLAT, grb$GLON)
ss
euc = euclid(ss)
euc
# Plot a wireframe sphere (unit sphere)
spheres3d(0, 0, 0, radius = 1, color = "gray", alpha = 0.1, front = "lines", back = "lines")

# Plot the data points, colored by cluster
colors <- c("red", "blue", "green")[clusters]  # Assign colors to clusters
points3d(euc[, 1], euc[, 2], euc[, 3], color = colors, size = 5)


# Add a legend
legend3d("topright", legend = paste("Cluster", 1:3), col = c("red", "blue", "green"), pch = 16)

# Add axes for reference
axes3d()

rgl.snapshot("my3d_skmeans.png")




cluster_2 = result$clusters
spheres3d(0, 0, 0, radius = 1, color = "gray", alpha = 0.1, front = "lines", back = "lines")

# Plot the data points, colored by cluster
colors <- c("red", "blue", "green")[cluster_2]  # Assign colors to clusters
points3d(euc[, 1], euc[, 2], euc[, 3], color = colors, size = 5)


# Add a legend
legend3d("topright", legend = paste("Cluster", 1:3), col = c("red", "blue", "green"), pch = 16)

# Add axes for reference
axes3d()

rgl.snapshot("my3d_skmeans.png")
